Local ecological divergence of two closely related stag beetles based on genetic, morphological, and environmental analyses

Abstract The process of phenotypic adaptation to the environments is widely recognized. However, comprehensive studies integrating phylogenetic, phenotypic, and ecological approaches to assess this process are scarce. Our study aims to assess whether local adaptation may explain intraspecific differentiation by quantifying multidimensional differences among populations in closely related lucanid species, Platycerus delicatulus and Platycerus kawadai, which are endemic saproxylic beetles in Japan. First, we determined intraspecific analysis units based on nuclear and mitochondrial gene analyses of Platycerus delicatulus and Platycerus kawadai under sympatric and allopatric conditions. Then, we compared differences in morphology and environmental niche between populations (analysis units) within species. We examined the relationship between morphology and environmental niche via geographic distance. P. kawadai was subdivided into the “No introgression” and “Introgression” populations based on mitochondrial COI gene – nuclear ITS region discordance. P. delicatulus was subdivided into “Allopatric” and “Sympatric” populations. Body length differed significantly among the populations of each species. For P. delicatulus, character displacement was suggested. For P. kawadai, the morphological difference was likely caused by geographic distance or genetic divergence rather than environmental differences. The finding showed that the observed mitochondrial–nuclear discordance is likely due to historical mitochondrial introgression following a range of expansion. Our results show that morphological variation among populations of P. delicatulus and P. kawadai reflects an ecological adaptation process based on interspecific interactions, geographic distance, or genetic divergence. Our results will deepen understanding of ecological specialization processes across the distribution and adaptation of species in natural systems.


| INTRODUC TI ON
How and why the diversity of life on earth increased over time are key research questions in ecology and biogeography (Blanquart et al., 2013;Cox et al., 2016;Futuyma & Antonovics, 1992;Savolainen et al., 2013;Thomas et al., 2016). Genetic and ecological speciation can occur in different parts of an ancestral species' range in which contrasting environmental conditions lead directly or indirectly to the evolution of reproductive isolation (Faulkes et al., 2004;Rundle & Nosil, 2005;Schluter, 2001). However, genetic divergence within and among species does not always cause divergence of morphological and other phenotypic traits due to silent mutations and phenotypic convergence (Fitch, 1970;Ujvari et al., 2015). Adaptative phenotypic variation often occurs via the evolution of eco-morphological innovations that contribute to ecological specialization in response to environmental variations or interspecific interactions (Devictor et al., 2010;Mammola et al., 2020). Therefore, evaluation of the phylogenetic constraints on traits and trait-environment relationships can elucidate the mechanisms underlying evolutionary selection and their impact on current ecological patterns.

Considering adaptation via multivariate genetic and trait analyses
is essential in such situations. However, comprehensive studies integrating phylogenetic, phenotypic, and ecological approaches to assessing speciation process and identifying phenotypic variations correlated with local adaptation have usually been neglected.
Here, we investigated inter-and intraspecific relationships using genetic, morphological, and ecological data for two closely related Platycerus beetles, Platycerus delicatulus Lewis, 1883, andPlatycerus kawadai Fujita andIchikawa, 1982, to explore how local adaptation shapes their habitat preference. P. delicatulus and P. kawadai of the family Lucanidae are endemic to Japan and exhibit geographic genetic variations (Kubota et al., 2011). Both species prefer mature cool temperate deciduous broad-leaved forests. P. delicatulus has a wide distribution across the main islands of Japan, except Hokkaido.
P. kawadai appears to be endemic to central Japan ( Figure 1). Both species co-occur throughout the range of P. kawadai, although some differences in host wood preference have been observed. More specifically, P. delicatulus and P. kawadai prefer hard and dry decaying wood as their larval resources, whereas all other Platycerus species in Japan prefer soft and wet decaying wood on the forest floor. However, P. delicatulus is more abundant at lower elevations, especially on thick decaying wood, and P. kawadai tends to target thin decaying wood at higher elevations . Two species would lose large portions of present suitable area under climate change (Zhang & Kubota, 2021). Phylogenetically, the two species diverged approximately 1 million years ago, and no hybridization between them has been recorded (Kubota et al., 2011;Zhu et al., 2020). P. delicatulus and P. kawadai are sister species with similar morphological and ecological attributes, such that sympatric distributions might lead to ecological divergence. Congeneric and ecologically similar species are considered good models for studies of ecological divergence, and thus these two species provide an opportunity to explore mechanisms underlying niche evolution and develop policies for insect management and conservation strategies.
The present study aimed to quantify multidimensional differences among populations that may arise due to local adaptation in the closely related species P. delicatulus and P. kawadai. First, we estimated the intra-and interspecific evolutionary dynamics of these two species across their ranges and constructed intraspecific analysis units using integrated phylogenetic results for F I G U R E 1 Occurrence records of Platycerus delicatulus (a) and Platycerus kawadai (b) at the collection sites in Japan both species under sympatric or allopatric conditions. We then explored differences in morphology and environmental niche among the populations within each species. We examined the relationship between morphology and environmental niche via geographic distance to assess whether local adaptation may explain population differentiation.

| Molecular procedures and phylogenetic analyses
This study focused on P. delicatulus and P. kawadai individuals collected from 2005 to 2020 for genetic analysis across the entire geographic range of these two species (Figure 1). The collection sites of the two species are listed in Appendix 1. Besides, Platycerus akitaorum Imura, 2007, andPlatycerus sugitai Okuda &Fujita, 1987, were used as outgroups.
In this study, we determined 94 and 45 sequences of the mitochondrial cytochrome oxidase subunit I (COI) gene and nuclear internal transcribed spacer (ITS) region, respectively (Appendix 2).
Genomic DNA was extracted from the testis or muscle tissues of adult beetles or larvae preserved in absolute ethyl alcohol using the Wizard Genomic DNA Purification kit (Promega).
We amplified fragments of the COI gene (primers C1-J-2183 and L2-N-3014, Simon et al., 1994) and ITS region (primers 5.8S38F and ITS4col, Tanahashi & Hawes, 2016) to explore the phylogenetic relationships within and between the two species. COI was amplified by polymerase chain reaction (PCR) at 94°C for 3 min, followed by 30 cycles of 94°C for 1 min, 48°C for 1 min, and 72°C for 1 min, and a final extension for 7 min at 72°C. The ITS region was amplified using the same process, but with an annealing temperature of 50°C. The PCR products were purified using the Illustra ExoStar Clean-Up kit (GE Healthcare).
Additionally, we used 65 COI and 5 ITS sequences for P. delicatulus and P. kawadai, and 9 COI and 2 ITS sequences for the outgroup (P. akitaorum and P. sugitai) from previous studies (Kubota et al., 2010(Kubota et al., , 2011Zhu et al., 2020). In total, we used 168 COI and 52 ITS sequences for analysis. The best-fit substitution model for COI and the ITS region were selected using jModelTest v.2.1.10 (Darriba et al., 2012) based on the Akaike information criterion (AIC).

F I G U R E 2
The eight investigated morphological traits investigated in this study. All traits were measured on the right side of the beetle's body, with the left side measured only when body parts were damaged or missing F I G U R E 3 Male genital endophallus of Platycerus delicatulus (a, c) and P. kawadai (b, d (Papadopoulou et al., 2010).
The data consisted of only in-group taxa, and the topology was fixed to the ML tree. Markov Chain Monte Carlo analysis was performed using 10 million generations, sampling every 1000 generations. The convergence of the chains was confirmed using Tracer v.

| Morphological analysis
For the morphological analysis, we assessed morphological external differentiation of P. delicatulus (central-to-northern Honshu) and performed a principal components analysis (PCA) using the procomp function in R v.3.6.3 (R Core Team, 2013) and visualized the results using "ggplot2" (Wickham & Wickham, 2007). To examine whether the two species differed in their morphological traits, we compared the principal component (PC) 1 and PC2 between two populations for each sex of each species. Then, we tested for body length (BL) differences between and within species and between the sexes using analysis of variance (ANOVA) and Tukey's HSD post hoc tests; BL was used as the measure for analysis, as it provides greater reproducibility than an axis derived from PCA (Barton et al., 2011).
For genital morphology, although we observed little difference in endophallic structure between P. delicatulus and P. kawadai (Figure 3), which may be concerning for reproductive isolation, we found no difference among populations within each species. Quantitatively assessing the membranous part of the endophallus is difficult, so we did not consider genital morphological variation.

| Environmental analysis
Environmental data were downloaded from the Worldclim database To quantify the environmental niches of P. delicatulus and P.
kawadai populations, we used two statistical approaches. First, PCA was performed on the environmental variables using procomp function in R v.3.6.3 (R Core Team, 2013) and visualized using "ggplot2" (Wickham & Wickham, 2007). Second, we compared the environmental niche spaces of the species using n-dimensional hypervolumes analyses (Hutchinson, 1957), which were conducted using the "hypervolume" R package (Blonder et al., 2018). We constructed the hypervolumes using the six retained variables for the major populations. All environmental variables were natural log-transformed for analysis. All hypervolumes were created using the Gaussian kernel density estimator method with the default Silverman bandwidth estimator (Blonder et al., 2014(Blonder et al., , 2018. To compare hypervolumes among environmental variables, we quantified the pairwise overlap between populations, using the Jaccard and Sorensen similarity indexes following Blonder et al. (2018).

| Correlations between morphology and environmental niche
We conducted Mantel tests and partial Mantel tests using the "vegan" R package to test correlation between the morphological and environmental distances of P. delicatulus and P. kawadai (Oksanen et al., 2013). Morphological distance was calculated as the Euclidean pairwise distance of BL between localities because BL is considered as an important trait for resource competition and reproductive interference (Okuzaki, 2021;

| Phylogenetic relationship between species
We sequenced 784 bp of the COI gene and 730-732 bp of the ITS region. These sequences were deposited in GenBank (DDBJ accession numbers: LC651809-LC651901 for the COI gene, and LC651902-LC651946 for the ITS region). The best-fit models were GTR + I + G for COI and GTR + G for the ITS region. including Site 97 exhibited no introgression type. Sites at which no genetic samples were collected were assigned to the category of the closest site at which genetic samples were collected. We subdivided P. delicatulus populations belonging to Clade II-a-2 into "Sympatric" population and "Allopatric" population. Sympatric population range covers whole range of P. kawadai, whereas both species cannot be always collected at the same site ( Figure 1, Appendix 1). In the following part, we examined the morphological differentiation among these analysis units of two species.

F I G U R E 9
Morphological differentiation between populations with respect to variations in body length (BL) for both male (a) and female (b) individuals. Analysis of variance (ANOVA) results are also shown. Different letters indicate significant differences between populations (Tukey's test: p < .05) For P. kawadai, PC1 explained 70.79% and 55.94% of the total variance for male and female, respectively. The significant difference between the populations in PC2 was detected only for P. delicatulus males (Figure 8). In this case, the eigenvalue of PC2 was 0.71 and the highest loading score for PC2 was 0.68 of head length (HL) ( Table 3).  (Figure 9a). Table 1 On the other hand, female BL varied significantly between the two species ( Figure 9b). Sympatric population of P. delicatulus and No introgression population of P. kawadai showed the highest and lowest value, respectively (Figure 9).

| Environmental niche
For P. delicatulus, we found the PCA results suggested that Sympatric population had a narrower environmental space than that of Allopatric population, especially in terms of elevation, temperature seasonality (Bio4), and mean temperature in wettest quarter (Bio8) ( Figure 10a;  Figure 10b). No introgression population exhibited higher temperature seasonality and lower mean temperature of wettest quarter, favoring less precipitation (Bio12 and Bio19) and a wider elevation compared with the Introgression population of P. kawadai (Table 5).
The multidimensional variations in the environmental space of both species are shown as niche hypervolumes in Figures 11 and 12, illustrating that the populations occupied different ecological spaces with relatively little overlap. For P. delicatulus, the niche hypervolume was much greater for the Allopatirc population than for the Sympatric population, and they overlapped slightly (Sørensen similarity = 0.057, Jaccard similarity = 0.029; Figure 11). For P. kawadai, the Sørensen and Jaccard similarity index values of the hypervolumes were 0.135 and 0.072 in No introgression and Introgression populations, respectively. Generally, Bio4 did not overlapped between the populations (Figure 12).

| Correlation between morphological and environmental niche
For P. delicatulus, simple Mantel tests showed that the morphological distance between populations was not significantly correlated with environmental (male, p = .104; female, p = .283) or geographic distances (male, p = .119; female, p = .315) ( Table 6). Morphological distance was not related with environmental distance after controlling for the effect of geographic distance (male, p = .102; female, p = .241, Figure 13a,c) and with geographic distance after controlling for environmental distance (male, p = .608; female, p = .588, Figure 13b,d)

F I G U R E 11
Hypervolumes obtained from multidimensional kernel density estimation of the studied population (Allopatric and Sympatric population) of Platycerus delicatulus based on weakly correlated environmental variables. The larger colored dots represent species centroids F I G U R E 1 2 Hypervolumes obtained from multidimensional kernel density estimation of the studied population (Allopatric and Sympatric population) of Platycerus kawadai based on weakly correlated environmental variables. The larger colored dots represent species centroids TA B L E 6 Single and partial Mantel test results based on morphological, environmental, and geographic distances between occurrence sites of Platycerus delicatulus and P. kawadai based on the partial Mantel test results. On the other hand, for P.
kawadai, morphological distance was significantly correlated with the environmental (male, p = .005; female, p = .03) and geographic distances (male, p < .001; female, p = .003) ( Table 6). Morphological distances were not significantly correlated with environmental distances after controlling for geographic distances using partial Mantel tests for P. kawadai (male, p = .470; female, p = .698, Figure 14a,c), however, morphological distance was significantly correlated with geographic distances after controlling for environmental distance in the same manner (male, p < .001; female, p = .010, Figure 14b,d).

| Phylogeographic history of the two related species
The genetic sample collection sites of the two species cover almost their entire distribution ranges (Appendix 2). Phylogenetic analyses based on the ITS region suggested that P. delicatulus and P. kawadai are each essentially monophyletic (Figure 4). This result aligns with the phylogenetic results of their yeast symbionts . Since the ancestral branches of P. delicatulus diverged in western Japan, it is likely that the two species were separated and speciated in western (P. delicatulus) and central (P. kawadai) Japan approximately 1.16 Mya (Figures 5 and 6).
Following that speciation event, P. delicatulus was separated into Based on our observation, males of Platycerus species always try F I G U R E 1 3 Partial regression plots illustrating the relationship between morphological distance and the environmental distance controlling geographic distance (a and c), and between morphological distance and geographic distance controlling for environmental distance (b and d) for male and female of Platycerus delicatulus, respectively to mate immediately with any female during the reproductive season. When there is a chance of heterospecific mating, interspecific differences in body size and genitalia size may work as premating and mechanical isolation mechanisms, respectively (Kubota & Sota, 1998;Okuzaki, 2021). A similar phylogeographic pattern has been documented in other beetles (Kosuda et al., 2016;Zhang & Sota, 2007).

These results indicated that No introgression population and
Introgression population of P. kawadai differed mainly in terms of COI, but they cannot be distinguished using ITS sequences.
Possible explanations for the mitochondrial-nuclear discordance could be associated with sex-biased dispersal, mating, and offspring production (Bonnet et al., 2017). Genetic drift is ubiquitous in populations and can interact with many of the above processes to increase discordance between mitochondrial and nuclear genes (Toews & Brelsford, 2012). But it is difficult to explain the essential topological difference between the COI and ITS phylogenies just for these reasons. Another possible evolutionary scenario for such a discordance is the incomplete lineage sorting following the ancestral polymorphism of mitochondrial gene (Funk & Omland, 2003). However, it is unlikely that the ancestor of P. kawadai had possessed both mitochondrial Clades I and II-a-1 because Clade IIa-1 had occurred in a P. delicatulus type subclade (Clade II-a) after initial geographical differentiation within P. delicatulus. An alternative and more likely scenario is historical mitochondrial introgression following the range expansion of these species. Because Clade II-a-1 was diverged from a P. delicatulus type clade around 0.74 Mya, the replacement by an introgressive clade seems to be very rare and only one replacement is recognized.

| Factors affecting morphological differences among Platycerus populations within species
In this study, we constructed intraspecific analysis units of two Platycerus species based on interspecific ranges and evolutionary dynamics, and then evaluated the factors affecting the morphological differences within each species. Among the eight morphological traits shown in Figure 3, BL was the most effective variable for explaining morphological variation (Figure 7). Meanwhile, the results of the n-dimensional hypervolume analysis revealed environmental heterogeneity among populations. We tested whether the morphological F I G U R E 1 4 Partial regression plots illustrating the relationship between morphological distance and the environmental distance controlling geographic distance (a and c), and between morphological distance and geographic distance controlling for environmental distance (b and d) for male and female of Platycerus kawadai, respectively. Lines represent significant regressions of the residuals variation across populations was better explained by geographic distance with dispersal or by environmental filtering for studied species.
For P. delicatulus, the morphological (BL) distance among collection sites was not correlated with environmental factors or with geographic distance, and therefore these factors could not explain the morphological divergence between Allopatric and Sympatric populations. The latter population is larger than the former, and likely arose via character displacement against P. kawadai ( Figure 9).
As P. delicatulus and P. kawadai are capable of mating, the putative character displacement may be caused by reproductive interference other than the resource competition. Overall, our results suggest that interspecific interaction has played a major role in driving the morphological differentiation of P. delicatulus populations.
For P. kawadai, morphological distance was correlated with geographic distance after controlling for environmental distance (Table 6). This result suggests that geographic distance (i.e., low dispersal ability) might have led to morphological differentiation.
Therefore, dispersal is assumed to drive the morphological diversification of populations. Meanwhile, dispersal ability could influence range limits and gene flow among populations, which may be associated with niche differentiation. In addition, previous studies showed that morphological adaptation to local ecology can also have resulted from phenotypic plasticity or from genetic differences among populations (Borokini et al., 2021;Ghalambor et al., 2007;Kunz et al., 2022;Price et al., 2003;Schmid & Guillaume, 2017). Although phenotypic plasticity has been documented in response to variations in multiple environmental variables (Chevin & Lande, 2015;Gratani, 2014;Lande, 2009;Wang et al., 2021), we found morphological distance was not correlated with environmental distance after controlling for geographic distance (Table 6). Thus, environmental factors are unlikely to be responsible for the observed morphological differentiation in P. kawadai. However, we cannot exclude the possibility that genetic divergence, such as that achieved via genetic drift and intra-and interspecific gene flow, promoted the morphological divergence. Further studies are required to verify whether this possibility would explain the morphological differentiation among populations of P. kawadai.
Populations often experience different environmental conditions, leading to the evolution of different phenotypes to maximize fitness (Freudiger et al., 2021;Jones et al., 2021). Most studies have shown that body size is affected by environmental filtering and food availability, which exhibit trade-off relationships (Dmitriew, 2011;Konuma et al., 2011;Runemark et al., 2015). Our results showed that intraspecific morphological variations in P. delicatulus and P. kawadai were related to interspecific interaction and geographic distance, respectively. These results indicated divergence between populations in directions of morphological variation and provided significant insights into species adaptation processes.
In conclusion, we integrated morphological, environmental, and molecular data across the geographic ranges of two species to investigate the ecological-evolutionary processes that may drive divergence processes among populations and across geography.
We found that morphological and ecological niche differentiation within species may be driven by interspecific interaction, as well as dispersal ability. These differentiations may associate with specialization for habitat preference. Our results elucidate ecological process across species' distributions through adaptation and plasticity in natural systems. Evidence of divergence between populations provides a useful reference for conservation strategies to enhance potential for adaptive response to the challenging climate changes.

ACK N OWLED G EM ENTS
We are grateful to N. Kubota, H. Otobe, S. Kobayashi, M. Ono, H.
Ono, E. Sakamoto, and M. Tanahashi for sample collection. This study was supported by the Japan Society for the Promotion of Science KAKENHI, Grant Numbers 25292082 and 21H04736 to K. Kubota.

CO N FLI C T O F I NTE R E S T S
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Sequence data are available at the DDBJ database under accession numbers LC651809-LC651901 for the COI gene, and LC651902-LC651946 for the ITS region, https://www.ddbj.nig.